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The exact solution of the BCS pairing Hamiltonian was found by 
Richardson in 1963. While little attention was paid to this exactly solv- 
able model in the remainder of the 20th century, there was a burst of 
work at the beginning of this century focusing on its applications in 
different areas of quantum physics. We review the history of this ex- 
act solution and discuss recent developments related to the Richardson- 
Gaudin class of integrable models, focussing on the role of these various 
models in nuclear physics. 



1. Cooper pairs, BCS and the Richardson exact solution 

The first breakthrough towards a microscopic description of the supercon- 
ducting phenomenon was due to CoopeiH, who in 1956 showed that a single 
pair of electrons on top of an inert Fermi sea could be bound by an infinites- 
imal attractive interaction. The search for a many-body wave function de- 
scribing a fraction of correlated and overlapping pairs mixed with a Fermi 
sea was a key goal for the rest of that year. Schrieffer came up with a 
solution at the beginning of 1957 and the BCS team (Bardeen, Cooper and 
Schrieffer) started an intensive and fruitful collaboration to explain quan- 
titatively many superconducting properties from the associated BCS wave 
function. This led to the famous BCS papeJ^ which provided a complete 
microscopic explanation of superconductivity. 

The success of the BCS theory quickly spread to other quantum many- 
body systems, including the atomic nucleus. In the summer of 1957, David 
Pines visited the Niels Bohr Institute and gave a series of seminars about 
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the yet unpublished BCS theory. Soon thereafter Bohr, Mottelson and 
Pines pubhshed a papei^^ suggesting that the gaps observed in even-even 
nuclei could be due to superconducting correlations. They noted, how- 
ever, that these effects should be strongly influenced by the finite size of 
the nucleus. Since then, and up to the present, number projection and 
in general symmetry restoration in the BCS and Hartree-Fock-Bogoliubov 
approximations have been important issues in nuclear structure. 

At the beginning of the sixties, while several groups were developing nu- 
merical techniques for number-projected BCS calculation^^, Richardson 
provided an exact solution for the reduced BCS HamiltoniadSE. In spite of 
the importance of his exact solution, this work did not have much impact 
in nuclear physics with just a few exceptions. Later on, his exact solution 
was rediscovered in the framework of ultrasmall superconducting grain^ 
where BCS and number-projected BCS were unable to describe appropri- 
ately the crossover from superconductivity to a normal metal as a function 
of the grain size. Since then, there has been a flurry of work extending the 
Richardson exact solution to families of exactly-solvable models, now called 
the Richardson-Gaudin (RG) model^, and applying these models to dif- 
ferent areas of quantum many-body physics including mesoscopic systems, 
condensed matter, quantum optics, cold atomic gases, quantum dots and 
nuclear structur^^^. In this paper, we review Richardson's solution, its gen- 
eralization to the exactly-solvable RG models and discuss the applications 
of these models in nuclear physics. 

2. The Richardson solution of the reduced BCS hamiltonian 

We will focus on a pairing Hamiltonian with constant strength G acting in 
a space of doubly-degenerate time- reversed states (k, fc). 

Hp = Y,('kclck - G^4clc^'Cfc' , (1) 

k k,k' 

where Ck are the single-particle energies for the doubly-degenerate orbits 

Cooper considered the addition of a pair of fermions with an attractive 
pairing interaction on top of an inert Fermi sea (FS) under the influence of 
this Hamiltonian. He showed that the pair eigenstate is 

l^Cooper) = E 44 l^'^) ' 
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where E is the energy eigenvalue. Cooper found that for any attractive 
value of G, the Fermi sea is unstable against the formation of such bound 
pairs. Therefore, an approach that takes into account a fraction of these 
correlated pairs mixed with a Fermi sea should be able to describe the 
superconducting phenomenon. 

The BCS approach followed a somewhat different path to the one sug- 
gested by Cooper, defining instead a variational wave function as a coherent 
state of pairs that are averaged over the whole system, 



l^scs) = e^' |0) , (3) 

where — J2k ^kclc^ is the coherent pair. Though errors due to the 
non-conservation of particle number in ([3]) are negligible when the number 
of pairs is sufficiently large, they can be important in such finite systems 
as atomic nucleF'. To accommodate these effects, number-projected BCS 
(PECS)'* considers a condensate of pairs of the form 

l^PBcs) = (r^'' |0) , (4) 

where M is the number of pairs and F^ has the same form as in BCS. 

Richardson^ proposed an ansatz for the exact solution of the pairing 
Hamiltonian ([T]) that followed closely Cooper's original idea. For a system 
with 2M + V particles, with v of these particles unpaired, his ansatz involves 
a state of the form 

\^)=bIbI--BIW) , (5) 

where the collective pair operators B^, have the form found by Cooper for 
the one-pair problem. 

Here L is the number of single-particle levels and 

W) = Wi,'^2 ■■■ (7) 

is a state of v unpaired fermions {v — i^^, with — \ or 0) defined by 
CkC-^\v) = 0, and nu \v) = Vk \v). 

In the one-pair problem, the quantities that enter ^ are the eigen- 
values of the pairing Hamiltonian, i. e. , the pair energies. Richardson pro- 
posed to use the M pair energies Ea in the many-body wave function of 
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Eqs. ([Slin]) as parameters which are chosen to fulfill the eigenvalue equation 
Hp 1 4') — E\'^). He showed that this is the case if the pair energies satisfy 
a set of M non-linear coupled equations 



which are now called the Richardson equations. The second term repre- 
sents the interaction between particles in a given pair and the third term 
represents the interaction between pairs. The associated eigenvalues of H 
are given by 



namely as a sum of the pair energies. 

Each independent solution of the set of Richardson equations defines a 
set of M pair energies that completely characterizes a particular eigenstate 
(l5l[6]). The complete set of eigenstates of the pairing Hamiltonian can be 
obtained in this way. The ground state solution is the energetically lowest 
solution in the v = oi v = 1 sector, depending on whether the system has 
an even or an odd number of particles, respectively. 

There are a couple of points that should be noted here. First, in contrast 
to the BCS solution, each Cooper pair i?J is distinct. Second, if one of 
the pair energies Ea is complex, then its complex-conjugate is also a 
solution. From this latter point we see that j^*) preserves time-reversal 
invariance. 

On inspection of the Richardson pair dH), we see that a pair energy that 
is close to a particular 2efc, i.e. close to the energy of an unperturbed pair, is 
dominated by this particular configuration and thus defines an uncorrelated 
pair. In contrast, a pair energy that lies sufficiently far away in the complex 
plane produces a correlated Cooper pair. This is to be contrasted with the 
single BCS coherent pair, which has amplitude Zk = Vk/uk and which mixes 
correlated and uncorrelated pairs over the whole system. 

3. Generalization to the Richardson-Gaudin class of inte- 
grable models 





L 



M 




(9) 



In this section, we discuss how to generalize the standard pairing model, 
which as we have seen is exactly solvable, to a wider variety of exactly- 
solvable models, the so-called Richardson-Gaudin modelJ^^^, all of which are 
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based on the SU{2) algebra. We first introduce the generators of SU{2), 
using a basis more familiar to nuclear structure, 

^" = M E - ) ' = E Kr = (K+y. (10) 



Here creates a fermion in single-particle state jm, jm denotes the 
time reverse of jm, and ftj = j + i is the pair degeneracy of orbit j. 
These operators fulfill the SU(2) algebra [K+,KJ,] = 2djj>K°, [K°,Kp] = 

We now consider a general set of L Hermitian and number-conserving 
operators that can be built up from the generators of SU{2) with linear 
and quadratic terms, 



(11) 



Following Gaudirf^, we then look for the conditions that the matrices 
X and Y must satisfy in order that the R operators commute with one 
another. It turns out that there are essentially two families of solutions, 
referred to as the rational and hyperbolic families, respectively. 



The rational family 



The hyberbolic family 



Xij — Yij — (12) 



X., = 2^ , r., = (13) 
m - Vj m - Vj 

Here the set of L parameters rji are free real numbers. 

The traditional pairing model is an example of the rational family. It 
can be obtained as a linear combination of the integrals of motion. Hp = 
Ej^j-Ri(ej): with r]j ^Ej. 

The complete set of eigenstates of the rational integrals of motion is 
given by the Richardson ansatz ((SJ [S]) ■ This fact led GaudirP^I to try to 
relate his integrable models to the BCS Hamiltonian without success. The 
proof of integrability of the BCS Hamiltonian was found later in ref.li^. We 
will not present the general solution of the two integrable families here, 
referring the reader to refs.'^^li'. 
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The key point is that any Hamihonian that can be expressed as a hnear 
combination of the R operators can be treated exactly using this method. 
In the foUowing sections, we discuss nuclear applications of the standard 
pairing model and of a new model based on the hyberbolic family. 

4. Applications of the Richardson solution to pairing in nu- 
clear physics 

Richardson himself started to explore analytically the exact solution in 
nuclear structure for few pairs outside a doubly-magic coi^^^. He also 
proposed a numerical method to solve the equations for systems with 
equidistant levelJ^^, a model that was subsequently used as a benchmark 
to test many-body approximationJi^. However, the first application of the 
Richardson solution to a real nuclear system was reported by Andersson 
and Krumlinde^^ in 1977. They studied the properties of high-spin states 
in ^^"^Dy using an oblate deformed oscillator potential and including the 
effects of pairing at several different levels of approximation. They com- 
pared the results when pairing was treated with the traditional BCS ap- 
proximation, when it was treated in PBCS approximation (using the saddle 
point approximation) and when it was treated exactly using the Richardson 
method. 

Following that early work, there were sporadic references to the Richard- 
son method but no realistic studies of atomic nuclei until just a few years 
ago. In 2007, Dussel et al^^ reported a systematic study of pairing corre- 
lations in the even Sm isotopes, from ^'^^Sm through ^^^Sm, using the self- 
consistent deformed Hartree Fock-|-BCS method. The calculations made 
use of the density-dependent Skyrme force, SLy4, and treated pairing cor- 
relations using a pairing force with constant strength G assuming axial 
symmetry and taking into account 11 major shells. 

Using the results at self-consistency to define the HF mean field, pairing 
effects within that mean field were then considered using the alternative 
number-conserving PBCS approach and the exact Richardson approach. 
In this way it was possible to directly compare the three approaches to 
pairing with the same pairing Hamiltonian, a primary focus of the study. 
It should be noted here that the Hilbert space dimensions associated with 
the residual neutron pairing Hamiltonian is of the order of 3.9 x 10^'^ for 
^^'^Sm, whereas the exact Richardson approach requires the solution of a 
coupled set of 46 non-linear equations. 

In the one semi-magic nucleus ^^'^Sm that was studied, the principal 
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correlation effects arise when projection is included, taking the system from 
one that is normal at the level of BCS to one with substantial pairing 
correlations. Treating pairing exactly provides a further modest increase in 
pairing correlations of about 0.3 MeV. In non-semi-magic nuclei, the effect 
on the pairing correlation energy of the exact solution is significantly more 
pronounced. While there too number projection provides a substantial 
lowering of the energy, it now misses about 1 MeV of the exact correlation 
energy that derives from the Richardson solution. 
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Fig. 1. Square of the wave function of the most collective Cooper pairs in -'^■''^Sin (de- 
noted El, E2, E3, E4, and E5) and the pairing tensor (BCS) versus the single-particle 
level i in the Hartrec-Fock basis The results are presented for the physical value of the 
pairing strength, G = 0.106 MeV 

As noted earlier, the Richardson prescription gives rise to distinct 
Cooper pairs with distinct structure. This is illustrated in Figure 1 where 
we compare the square of the wave function for the most correlated Cooper 
pairs in ^^"^Sm, i.e. those whose pair energies lie farthest from the real 
axis in the complex plane, with the square of the pairing tensor UiVi that 
derives from the corresponding BCS solution. All wave functions are plot- 
ted versus the order of the single-particle states to make clear the relevant 
mixing of configurations in each pair. The pair label Ei refers to the two 
most collective pairs (with complex conjugate pair energies). E2 refers to 
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the next two most collective pairs, which are however only marginally col- 
lective. E3 refers to the next two in descending order of collectivity, but 
they only involve perturbative mixing of configurations and are not truly 
collective. The final two that are shown, E4 and E5, have real pair energies 
and involve almost pure single-particle configurations. ^From the figure, we 
see that even the most collective Cooper pairs are much less collective than 
UiVi, and therefore that their size in coordinate space is significantly larger 
than that of the BCS pairing tensoi'^, which is often used in the literature 
as a definition of the Cooper pair wave functiorPS. 

The exact Richardson solution was also used to study the gradual emer- 
gence of superconductivity in the Sn isotoped^. By making use of an exact 
mapping between the Richardson equations and a classical electrostatic 
problem in two dimensions, it was possible to get a physical picture of how 
superconductivity develops as a function of the pairing strength. In partic- 
ular, as the pairing strength is increased the pair energies gradually merge 
into larger structures in the complex plane as pair correlations gradually 
overcome single-particle effects. 

More recently, the Richardson solution has been applied to the treat- 
ment of pair correlations involving the continuum. The first work by 
Hasegawa and Kaneko-'^ considered only the effect of resonances in the 
continuum and as a result obtained complex energies even for the bound 
states of the system. Subsequent work by Id BetaiP^'^SI included the effects 
of the true continuum. The most recent pape^ treated nuclear chains that 
include both bound and unbound systems, e.g. the even-A Carbon isotopes 
up to ^^C. When the system is bound, the pair energies that contribute 
to the ground state occur in complex conjugate pairs, thus preserving the 
real nature of the ground state energy. Once the system becomes unbound 
this ceases to be the case. Now the pair energies that contribute to the 
ground state do not occur in complex conjugate pairs, explaining how a 
width arises in the energy of an unbound system within the Richardson 
approach. 

5. The hyperbolic model 

The hyperbolic family of models did not find a physical realization until 
very recently when it was shown that they could model a p-wave pairing 
Hamiltonian in a 2-dimensional lattic^^, such that it was possible to study 
with the exact solution an exotic phase diagram having a non-trivial topo- 
logical phase and a third-order quantum phase transition^^. Immediately 
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thereafter, it was shown that the hyperbohc family gives rise to a sepa- 
rable pairing Hamiltonian with 2 free parameters that can be adjusted to 
reproduce the properties of heavy nuclei as described by a Gogny HFB 
treatment!^. Both applications are based on a simple linear combination of 
hyperbolic integrals which give rise to the separable pairing Hamiltonian 

= ^ 77,K° - G ^ ^Tu^K+K^ . (14) 

If we interpret the parameters r]i as single-particle energies correspond- 
ing to a nuclear mean-field potential, the pairing interaction has the un- 
physical behavior of increasing in strength with energy. In order to reverse 
this unwanted effect, we define rji = 2{ei — a), where the free parameter 
a plays the role of an energy cutoff and Ei is the single-particle energy of 
the mean-field level i. Making use of the pair representation of SU{2), the 
exactly-solvable pairing Hamiltonian (fT4|) takes the form 

if = ^£,; (4c, + 4q) - 2G^V(a 44^1' c^', (15) 
i a' 

with eigenvalues E — 2aM + J2i^i'^i + "^p^P- The pair energies Ep 
correspond to a solution of the set of non-linear Richardson equations 

where Q = 2?J~T"'"^^~1- Each particular solution of Eq. ([T6| defines a 
unique eigenstate. 

Due to the separable character of the hyperbolic Hamiltonian, in BCS 
approximation the gaps = 2G^a — Si — Ei'Ui'Vi' — A^/a — Si 

and the pairing tensor UjVj = — A^a-e, have a very restricted 



form. In order to test the validity of the exactly solvable Hamiltonian (fT5|) 
we take the single-particle energies Ei from the HF energies of a Gogny HFB 
calculation and we fit the parameters a and G to the gaps and pairing ten- 
sor in the HF basis. Figure [5] shows the comparison for protons in 238jj 
between the Gogny HFB results in the HF basis and the BCS approxima- 
tion of the hyperbolic model. From these results we extracted the values 
a — 25.25 MeV and G = 2 x lO^'^MeV. The valence space determined 
by the cutoff a corresponds to 148 levels with 46 proton pairs. The size of 
the Hamiltonian in this space is 4.83 x 10'^*, well beyond the limits of exact 
diagonalization. However, the integrability of the hyperbolic model pro- 
vides an exact solution by solving a set of 46 non-linear coupled equations. 
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Moreover, the exact solution shows a gain in correlation of more than 2 
MeV suggesting the importance of taking into account correlations beyond 
mean-field. 




-50 -40 -30 -20 -10 ^0 10 20 



Fig. 2. Gaps A; and pairing tensor UiVi for protons in ^^^U. Open circles are Gogny 
HFB results in MeV. Solid lines are BCS results of the hyperbolic Hamiltonian in MeV. 



6. Extensions to non-compact and higher rank algebras 

Up to now, we have restricted our discussion to RG models that are based 
on the compact rank-1 SU{2) pair algebra. The method of constructing RG 
models can be extended to the non-compact rank-1 SU{1, 1) algebra as well, 
whereby pairing in bosonic systemd^ is described in complete analogy with 
the SU{2) case. An early application to the SO{Q) to [/(5) line of integra- 
bility of the Interacting Boson Model (IBM) was reported in ref.l^, with the 
exact solution being obtained there directly using an infinite dimensional 
algebraic technique. Further work on the IBM using the integrablc 5?7(1, 1) 
RG modef^ including high-spin bosons {d, g, ■ ■ ■) revealed a particular fea- 
ture of the repulsive boson pairing interaction that seems to provide a new 
mechanism for the enhancement of s — d dominance, giving further support 
for the validity of the s ~ d Interacting Boson Model. 

The RG models are not constrained to rank-1 algebras. They can be 
extended to any semi-simple Lie algebra^^. Richardson himself studied 
some restricted solutions of the T=l pairing modeP and the T=0,1 pair- 
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ing modeP^. As a general statement, the reduced pairing Hamiltonian is 
exactly solvable for any multi-component system. The first step in finding 
an exact solution is to identify the Lie algebra of the commuting pair op- 
erators and then to specialize the general solution given ir!22l. One has to 
keep in mind that while the SU (2) RG model has a single set of unknown 
parameters, the pair energies, larger rank algebras have as many sets of 
unknown parameters as the rank of the algebra. Therefore, the higher the 
rank of the algebra, the greater is the complexity of the solution. Several 
pairing Hamiltonians with relevance to nuclear physics have been studied 
in the last few years. 

i. The rank-2 5*0(5) RG modeP^ describes T=l proton-neutron pairing 
with non-degenerate single particle levels. The exact solution has two sets 
of spectral parameters, the pair energies and a second set associated with 
the SU (2) isospin subalgebra. In spite of the greater complexity, it was 
possible to solve exactly a T=l pairing Hamiltonian for the nucleus ^^Ge 
using a ^°Ca core, with a Hilbert space dimension well beyond the limits 
of exact diagonalization. 

ii. The rank-3 50(6) RG model'^^ describes color pairing, i. e. pairing 
between three-component fermions. The exact Richardson equations have 
three sets of spectral parameters, of which one correspond the the pair 
energies and the other two are responsible for the different couplings within 
the SU{3) color subalgebra. The model has been used to study the phase 
diagram of polarized three-component fermion atomic gases. However, it 
could in principle be exploited to describe non-relativistic quark systems. 

in. With increasing complexity, the rank-4 50(8) RG model^^ describes 
either T=0,1 proton- neutron pairing or four-component fermion gases. It 
contains four sets of spectral parameters. The model has been used to study 
alpha-like structures represented by clusters in the parameter space, and 
how these clusters dissolve into like-particle pairs with increasing isospin. 

iv. The rank-2 non-compact 5*0(3, 2) algebra generalizes the bosonic 
RG models to systems of interacting proton and neutron bosons^^. The 
model describes the IBM2 in the line of integrability between vibrational 
and 7-soft nuclei. The exact solution has been employed to study the 
influence of high-spin / and g bosons in the low-energy spectrum 

7. Summary and future outlook 

In this article, we have reviewed Richardson's solution of the pairing model 
and have discussed its generalization to a wider class of exactly solvable 
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models. Wc have also discussed the application of these models to a variety 
of problems in nuclear structure physics in which pairing plays a role. It 
should be noted here, however, that all of the models that we have discussed 
arc restricted to the pairing degree of freedom and thus do not allow explicit 
treatment of deformation effects. It is only through the use of Nilsson or 
deformed Hartree Fock single-particle energies that effects of deformation 
arc simulated. 

A key feature of the Richardson-Gaudin integrable models, is that they 
transform the diagonalization of the hamiltonian matrix, whose dimension 
grows exponentially with the size of the system, to the solution of a set 
of M coupled non-linear equations where Af is the number of pairs. This 
makes it possible to treat problems that could otherwise not be treated 
and in doing so to obtain information that is otherwise inaccessible. For 
example, we rciportcid an application of the rational RG pairing model to 
the even-mass Sm isotopes, where the size of the Hilbert space would exceed 
10^^ states, and an application of the hyperbolic RG pairing model to U^^^, 
where the size of the Hilbert space would exceed 10^^ states. In both cases, 
substantial gains in correlation energy were found when the problem was 
treated exactly. 

The exactly solvable RG Hamiltonians also provide excellent bench- 
marks for testing approximations beyond HFB in realistic situations both 

for even-even and odd-mass nuclei. Moreover, a self-consistent HF plus 
exact pairing approach could in principle be implemented to describe large 
regions of the table of nuclides. It might be possible to extend such a self- 
consistent approach to the 0(5) RG model, providing in this way a better 
description of those nuclei with N ^ Z in which T=l proton- neutron pair- 
ing correlations are expected to play a significant role. Unfortunately, the 
S0(8) T=0,1 RG model cannot accommodate the spin-orbit splitting in the 
single-particle energies. Nevertheless, this model could play an important 
role in helping to understand quartet clusterization and quartet condensa- 
tion in nuclear and cold atom systems. Finally, extension of the RG models 
to include the effects of the continuum seems to be an especially promising 
avenue to explore the physics of weakly-bound nuclei. 
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